In [ ]:
from __future__ import division, print_function
import os
import sys
from six.moves import cPickle as pickle
# Third-party
import astropy.coordinates as coord
import astropy.units as u
uno = u.dimensionless_unscaled
import matplotlib as mpl
import matplotlib.pyplot as pl
pl.style.use('apw-notebook')
%matplotlib inline
import numpy as np
# Custom
importgala.coordinates as gc
importgala.dynamics as gd
importgala.integrate as gi
importgala.potential as gp
fromgala.units import galactic
from scipy.signal import argrelmin
import ophiuchus.potential as op
from ophiuchus.data import OphiuchusData
from ophiuchus.util import integrate_forward_backward
from ophiuchus.coordinates import Ophiuchus
from ophiuchus import galactocentric_frame, vcirc, vlsr, RESULTSPATH
plotpath = "/Users/adrian/projects/ophiuchus-paper/figures/"
if not os.path.exists(plotpath):
os.mkdir(plotpath)
In [ ]:
ophdata = OphiuchusData()
ophdata_fit = OphiuchusData("(source == b'Sesar2015a') | (Name == b'cand9') | (Name == b'cand14')")
ophdata_fan = OphiuchusData("(source == b'Sesar2015b') & (Name != b'cand9') & (Name != b'cand14')")
all_names = ["static_mw"] + ["barred_mw_{}".format(i) for i in range(1,10)]
short_names = ["static"] + ["bar{}".format(i) for i in range(1,10)]
name_map = dict(zip(all_names, short_names))
Added 29 April 2016
In [ ]:
names = ['name', 'ra', 'dec', 'vlos', 'vlos_err']
all_bhb = np.genfromtxt("/Users/adrian/projects/ophiuchus/allstars.txt", usecols=range(5), names=names, dtype=None)
all_bhb_c = coord.ICRS(ra=all_bhb['ra']*u.degree, dec=all_bhb['dec']*u.degree)
all_bhb_c = all_bhb_c.transform_to(coord.Galactic)
In [ ]:
# global style stuff
orbit_style = dict(marker=None, color='k', alpha=0.05)
data_style = dict(marker='o', ms=4, ls='none', alpha=0.9, color='#2166AC',
markeredgecolor='k', markeredgewidth=0.5)
data_b_style = data_style.copy()
data_b_style['color'] = "#2166AC"
data_b_style['marker'] = "s"
fig,axes = pl.subplots(2,1,figsize=(4,4.5),sharex=True,sharey='row')
name = 'static_mw'
axes[1].set_xlabel("$l$ [deg]", fontsize=18)
path = os.path.join(RESULTSPATH, name, 'orbitfit')
w0 = np.load(os.path.join(path, 'w0.npy'))[:128].T
pot = op.load_potential(name)
orbit = integrate_forward_backward(pot, w0, t_forw=20., t_back=-20)
orbit_c,orbit_v = orbit.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame,
vcirc=vcirc, vlsr=vlsr)
orbit_l = orbit_c.l.wrap_at(180*u.deg)
orbit_oph = orbit_c.transform_to(Ophiuchus)
vr = (orbit_v[2].to(u.km/u.s)).value
# sky
axes[0].plot(ophdata_fit.coord.l.degree, ophdata_fit.coord.b.degree, **data_style)
axes[0].plot(ophdata_fan.coord.l.degree, ophdata_fan.coord.b.degree, **data_b_style)
axes[0].plot(all_bhb_c.l.degree, all_bhb_c.b.degree, ls='none', color='#666666', marker='o', alpha=0.4)
axes[0].yaxis.set_ticks(np.arange(27,32+1))
# radial velocity
axes[1].plot(ophdata_fit.coord.l.degree, ophdata_fit.veloc['vr'].to(u.km/u.s).value, **data_style)
axes[1].plot(ophdata_fan.coord.l.degree, ophdata_fan.veloc['vr'].to(u.km/u.s).value, **data_b_style)
axes[1].plot(all_bhb_c.l.degree, all_bhb['vlos'], ls='none', color='#666666', marker='o', alpha=0.4)
# axes[1].yaxis.set_ticks(np.arange(-300,300+1,100)) # 1
axes[1].yaxis.set_ticks(np.arange(225,325+1,25)) # 2
axes[0].set_xlim(9,2)
axes[0].set_ylabel("$b$ [deg]", fontsize=18)
axes[0].set_ylim(26.5, 32.5)
axes[1].set_ylabel(r"$v_r$ [${\rm km}\,{\rm s}^{-1}$]", fontsize=18)
# axes[1].set_ylim(-250, 350) # 1
axes[1].set_ylim(200, 350) # 2
# fig.tight_layout()
fig.subplots_adjust(left=0.3, right=0.98, top=0.96, bottom=0.15)
fig.savefig("/Users/adrian/projects/talks/thesis_colloquium/ophiuchus2.png", dpi=400)
In [ ]:
# global style stuff
orbit_style = dict(marker=None, color='#2166AC', alpha=0.05)
data_style = dict(marker='o', ms=6, ls='none', ecolor='#333333', alpha=0.75)
data_b_style = data_style.copy()
data_b_style['color'] = "#666666"
data_b_style['marker'] = "s"
fig,axes = pl.subplots(3,2,figsize=(6,7.5),sharex=True,sharey='row')
for i,name in enumerate(all_names[:2]):
axes[0,i].set_title(name_map[name], fontsize=20)
axes[2,i].set_xlabel("$l$ [deg]", fontsize=18)
axes[0,i].set_aspect('equal')
path = os.path.join(RESULTSPATH, name, 'orbitfit')
w0 = np.load(os.path.join(path, 'w0.npy'))[:128].T
pot = op.load_potential(name)
orbit = integrate_forward_backward(pot, w0, t_forw=20., t_back=-20)
orbit_c,orbit_v = orbit.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame,
vcirc=vcirc, vlsr=vlsr)
orbit_l = orbit_c.l.wrap_at(180*u.deg)
orbit_oph = orbit_c.transform_to(Ophiuchus)
vr = (orbit_v[2].to(u.km/u.s)).value
# sky
_tmp = data_style.copy(); _tmp.pop('ecolor')
axes[0,i].plot(ophdata_fit.coord.l.degree, ophdata_fit.coord.b.degree, **_tmp)
_tmp = data_b_style.copy(); _tmp.pop('ecolor')
axes[0,i].plot(ophdata_fan.coord.l.degree, ophdata_fan.coord.b.degree, **_tmp)
axes[0,i].plot(orbit_l.degree, orbit_c.b.degree, **orbit_style)
axes[0,i].yaxis.set_ticks(np.arange(27,32+1))
# distance
axes[1,i].errorbar(ophdata_fit.coord.l.degree, ophdata_fit.coord.distance.to(u.kpc).value,
ophdata_fit.coord_err['distance'].to(u.kpc).value, **data_style)
axes[1,i].errorbar(ophdata_fan.coord.l.degree, ophdata_fan.coord.distance.to(u.kpc).value,
ophdata_fan.coord_err['distance'].to(u.kpc).value, **data_b_style)
axes[1,i].plot(orbit_l.degree, orbit_c.distance.to(u.kpc).value, **orbit_style)
axes[1,i].yaxis.set_ticks(np.arange(6,9+1))
# radial velocity
axes[2,i].errorbar(ophdata_fit.coord.l.degree, ophdata_fit.veloc['vr'].to(u.km/u.s).value,
ophdata_fit.veloc_err['vr'].to(u.km/u.s).value, **data_style)
axes[2,i].errorbar(ophdata_fan.coord.l.degree, ophdata_fan.veloc['vr'].to(u.km/u.s).value,
ophdata_fan.veloc_err['vr'].to(u.km/u.s).value, **data_b_style)
axes[2,i].plot(orbit_l.degree, np.vstack(vr), **orbit_style)
axes[2,i].yaxis.set_ticks(np.arange(230,320+1,30))
axes[0,0].set_xlim(9,2)
axes[0,0].set_ylabel("$b$ [deg]", fontsize=18)
axes[0,0].set_ylim(26.5, 32.5)
axes[1,0].set_ylabel(r"$d_\odot$ [kpc]", fontsize=18)
axes[1,0].set_ylim(5.5, 9.5)
axes[2,0].set_ylabel(r"$v_r$ [${\rm km}\,{\rm s}^{-1}$]", fontsize=18)
axes[2,0].set_ylim(225, 325)
fig.tight_layout()
fig.savefig(os.path.join(plotpath, "orbitfits.pdf"))
fig.savefig(os.path.join(plotpath, "orbitfits.png"), dpi=400)
In [ ]:
# global style stuff
orbit_style = dict(marker=None, color='#2166AC', alpha=0.05)
data_style = dict(marker='o', ms=6, ls='none', ecolor='#333333', alpha=0.75)
data_b_style = data_style.copy()
data_b_style['color'] = "#666666"
data_b_style['marker'] = "s"
fig,axes = pl.subplots(2,2,figsize=(6,6),sharex=True,sharey='row')
for i,name in enumerate(all_names[:2]):
axes[0,i].set_title(name_map[name], fontsize=20)
axes[1,i].set_xlabel("$l$ [deg]", fontsize=18)
path = os.path.join(RESULTSPATH, name, 'orbitfit')
w0 = np.load(os.path.join(path, 'w0.npy'))[:128].T
pot = op.load_potential(name)
orbit = integrate_forward_backward(pot, w0, t_forw=20., t_back=-20)
orbit_c,orbit_v = orbit.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame,
vcirc=vcirc, vlsr=vlsr)
orbit_l = orbit_c.l.wrap_at(180*u.deg)
orbit_oph = orbit_c.transform_to(Ophiuchus)
mul = galactic.decompose(orbit_v[0]).value
mub = galactic.decompose(orbit_v[1]).value
# proper motion
axes[0,i].errorbar(ophdata_fit.coord.l.degree, galactic.decompose(ophdata_fit.veloc['mul']).value,
galactic.decompose(ophdata_fit.veloc_err['mul']).value, **data_style)
axes[0,i].plot(orbit_l.degree, np.vstack(mul), **orbit_style)
# axes[0,i].yaxis.set_ticks(np.arange(230,320+1,30))
# mub
axes[1,i].errorbar(ophdata_fit.coord.l.degree, galactic.decompose(ophdata_fit.veloc['mub']).value,
galactic.decompose(ophdata_fit.veloc_err['mub']).value, **data_style)
axes[1,i].plot(orbit_l.degree, np.vstack(mub), **orbit_style)
axes[0,0].set_xlim(9,2)
axes[0,0].set_ylim(-12,-2)
axes[1,0].set_ylim(-2,8)
axes[0,0].set_ylabel(r"$\mu_l$ [${\rm mas}\,{\rm yr}^{-1}$]", fontsize=18)
axes[1,0].set_ylabel(r"$\mu_b$ [${\rm mas}\,{\rm yr}^{-1}$]", fontsize=18)
fig.tight_layout()
fig.savefig(os.path.join(plotpath, "orbitfits-pm.pdf"))
fig.savefig(os.path.join(plotpath, "orbitfits-pm.png"), dpi=400)
In [ ]:
split_ix = 350
every = 50
mean_w0s = np.zeros((len(all_names), 6))
for i,name in enumerate(all_names):
with open(os.path.join(RESULTSPATH, name, 'orbitfit', 'sampler.pickle'), 'rb') as f:
sampler = pickle.load(f)
_x0 = np.vstack(sampler.chain[:,split_ix::every,:5])
mean_x0 = np.mean(_x0, axis=0)
std_x0 = np.std(_x0, axis=0)
transforms = [
lambda x: np.degrees(x),
lambda x: x,
lambda x: (x*u.radian/u.Myr).to(u.mas/u.yr).value,
lambda x: (x*u.radian/u.Myr).to(u.mas/u.yr).value,
lambda x: (x*u.kpc/u.Myr).to(u.km/u.s).value
]
cols = []
for j,_mean,_std in zip(range(len(mean_x0)), mean_x0, std_x0):
cols.append("{:.3f} {:.3f}".format(transforms[j](_mean), transforms[j](_std)))
print(" & ".join(cols))
mul = (mean_x0[2]*u.radian/u.Myr).to(u.mas/u.yr).value
mean_w0s[i] = ophdata._mcmc_sample_to_w0(mean_x0)[:,0]
In [ ]:
split_ix = 350
every = 50
for i,name in enumerate(all_names):
with open(os.path.join(RESULTSPATH, name, 'orbitfit', 'sampler.pickle'), 'rb') as f:
sampler = pickle.load(f)
_x0 = np.vstack(sampler.chain[:,split_ix::every,5:])
mean_x0 = np.mean(_x0, axis=0)
print("{:.2f} {:.2f} {:.2f}".format((mean_x0[0]*u.radian).to(u.deg), mean_x0[1], (mean_x0[2]*u.kpc/u.Myr).to(u.km/u.s)))
std_x0 = np.std(_x0, axis=0)
print("{:.2f} {:.2f} {:.2f}".format((std_x0[0]*u.radian).to(u.deg), std_x0[1], (std_x0[2]*u.kpc/u.Myr).to(u.km/u.s)))
print()
In [ ]:
_tmp_cache = dict()
In [ ]:
fig,axes = pl.subplots(2,5,figsize=(9,5),sharex=True,sharey=True)
for i,name in enumerate(all_names):
this_w0 = mean_w0s[i]
pot = op.load_potential(name)
if name not in _tmp_cache:
print("integrating")
orbit = pot.integrate_orbit(this_w0, dt=-1., nsteps=6000., Integrator=gi.DOPRI853Integrator)
_tmp_cache[name] = orbit
else:
orbit = _tmp_cache[name]
print(orbit.pericenter(), orbit.apocenter())
axes.flat[i].plot(orbit.pos[1], orbit.pos[2], marker=None)
axes.flat[i].set_title(name_map[name], fontsize=18)
if i > 4:
axes.flat[i].set_xlabel("$y$ [kpc]", fontsize=18)
axes[0,0].set_ylabel("$z$ [kpc]", fontsize=18)
axes[1,0].set_ylabel("$z$ [kpc]", fontsize=18)
_s = 17
axes[0,0].set_xlim(-_s,_s)
axes[0,0].set_ylim(-_s,_s)
axes[0,0].xaxis.set_ticks([-10,0,10])
axes[0,0].yaxis.set_ticks([-10,0,10])
fig.tight_layout()
fig.savefig(os.path.join(plotpath, "orbit-yz.png"), dpi=300)
fig.savefig(os.path.join(plotpath, "orbit-yz.pdf"))
In [ ]:
for i,name in enumerate(all_names):
orbit = _tmp_cache[name]
pl.figure()
pl.plot(orbit.t, np.sqrt(np.sum(orbit.pos**2,axis=0)))
pl.plot(orbit.t, np.abs(orbit.pos[2]))
pl.xlim(-600,10)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
# global style stuff
orbit_style = dict(marker=None, color='#2166AC', alpha=0.05)
data_style = dict(marker='o', ms=4, ls='none', ecolor='#666666', alpha=0.75)
for n,name_subset in enumerate([all_names[:5], all_names[5:]]):
fig,axes = pl.subplots(3,5,figsize=(9,6.5),sharex=True,sharey='row')
for i,name in enumerate(name_subset):
axes[0,i].set_title(name_map[name], fontsize=20)
axes[2,i].set_xlabel("$l$ [deg]", fontsize=18)
path = os.path.join(RESULTSPATH, name, 'orbitfit')
w0 = np.load(os.path.join(path, 'w0.npy'))[:128].T
pot = op.load_potential(name)
orbit = integrate_forward_backward(pot, w0, t_forw=16., t_back=-10)
orbit_c,orbit_v = orbit.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame,
vcirc=vcirc, vlsr=vlsr)
orbit_oph = orbit_c.transform_to(Ophiuchus)
vr = (orbit_v[2].to(u.km/u.s)).value
# sky
_tmp = data_style.copy()
_tmp.pop('ecolor')
axes[0,i].plot(ophdata.coord.l.degree, ophdata.coord.b.degree, **_tmp)
axes[0,i].plot(orbit_c.l.degree, orbit_c.b.degree, **orbit_style)
axes[0,i].yaxis.set_ticks(np.arange(27,32+1))
# distance
axes[1,i].errorbar(ophdata.coord.l.degree, ophdata.coord.distance.to(u.kpc).value,
ophdata.coord_err['distance'].to(u.kpc).value, **data_style)
axes[1,i].plot(orbit_c.l.degree, orbit_c.distance.to(u.kpc).value, **orbit_style)
axes[1,i].yaxis.set_ticks(np.arange(6,9+1))
# radial velocity
axes[2,i].errorbar(ophdata.coord.l.degree, ophdata.veloc['vr'].to(u.km/u.s).value,
ophdata.veloc_err['vr'].to(u.km/u.s).value, **data_style)
axes[2,i].plot(orbit_c.l.degree, np.vstack(vr), **orbit_style)
axes[2,i].yaxis.set_ticks(np.arange(230,320+1,30))
axes[0,0].set_xlim(9,2)
axes[0,0].set_ylabel("$b$ [deg]", fontsize=18)
axes[0,0].set_ylim(26.5, 32.5)
axes[1,0].set_ylabel(r"$d_\odot$ [kpc]", fontsize=18)
axes[1,0].set_ylim(5.5, 9.5)
axes[2,0].set_ylabel(r"$v_r$ [${\rm km}\,{\rm s}^{-1}$]", fontsize=18)
axes[2,0].set_ylim(225, 325)
fig.tight_layout()
# fig.savefig(os.path.join(plotpath, "orbitfits-{}.pdf".format(n)))
# fig.savefig(os.path.join(plotpath, "orbitfits-{}.png".format(n)))